home *** CD-ROM | disk | FTP | other *** search
/ Disc to the Future 2 / Disc to the Future Part II Programmer's Reference (Wayzata Technology)(6013)(1992).bin / MAC / THINKC / 4_0 / ORBIT_SO / ORBIT.C < prev    next >
C/C++ Source or Header  |  1992-07-19  |  5KB  |  250 lines
  1. /* Copyright 1987                            */
  2. /* David Palmer                             */
  3. /* Mail code 220-47                            */
  4. /* California Institute Of Technology         */
  5. /* Uses EasyDialog.c  (also ⌐ 1987 By David Palmer) */
  6. /* Duplication, modification, and examination allowed on a    */
  7. /* non-commercial basis only.  Commercial use prohibited    */
  8. /* without prior written agreement with the author.  (This    */
  9. /* includes sale by for-profit companies, and use as an        */
  10. /* inducement to buy something.)                            */
  11.  
  12. #include <quickdraw.h>
  13. #include <WindowMgr.h>
  14. #include "orbit.h"
  15. #include <EventMgr.h>
  16.  
  17. #include <stdio.h>
  18.  
  19. #define G 6.67e-11            /* MKS */
  20. #define NMAX 100            /* how many particles */
  21. #define TRUE (-1)
  22. #define FALSE 0
  23.  
  24. #define HUGE 1e4000
  25. double sqrt();
  26.  
  27. double precision = 100;
  28. double drawperiod = 86400;
  29. double scale = 2083333333;
  30. int trailstyle = 3;
  31. int fcenter = TRUE;
  32. int background = 1;
  33.  
  34. PARTICLE rgpa[NMAX];        /* All of the particles */
  35.  
  36. int cpa;                        /* number of particles */
  37. int cparead;
  38. int xzero;
  39. int yzero;
  40. int drawcode, erasecode;
  41.  
  42. double dt = 864;
  43. double rmin;
  44. double v2max;
  45.  
  46. main(argc, argv)
  47. int argc;
  48. char *argv[];
  49. {
  50.     InitGraf (&thePort);    /* initialize the screen         */
  51.     InitFonts();        
  52.     InitWindows();
  53.     InitMenus();
  54.     InitCursor();        /* make the arrow Cursor         */
  55.     TEInit();
  56.     InitDialogs(0L);
  57.  
  58.     DoIt();
  59. }
  60.  
  61. CenterParticles(ppa, cpa)
  62. PARTICLE *ppa;
  63. int cpa;
  64. {
  65.     int i, j;
  66.     double p[NDIMS];        /* total momentum    */
  67.     double mx[NDIMS];        /* first moment        */
  68.     double m=0;
  69.     
  70.     for (j = 0 ; j < NDIMS ; j++) {
  71.         p[j] = 0;
  72.         mx[j] = 0;
  73.     }
  74.     for (i = 0 ; i < cpa ; i++) {
  75.         m += ppa[i].m;
  76.         for (j = 0 ; j < NDIMS ; j++) {
  77.             mx[j] += ppa[i].m * ppa[i].x[j];
  78.             p[j] += ppa[i].p[j];
  79.         }
  80.     }
  81.     if (m == 0)            /* can happen with negative and positive masses */
  82.         return;
  83.     for (j = 0 ; j < NDIMS ; j++) {
  84.         mx[j] /= m;
  85.         p[j] /= m;                /* velocity of center of mass */
  86.     }
  87.     for (i = 0; i < cpa ; i++)
  88.         for (j = 0 ; j < NDIMS ; j++) {
  89.             ppa[i].x[j] -= mx[j];
  90.             ppa[i].p[j] -= ppa[i].m * p[j];
  91.             ppa[i].xt[j] = ppa[i].x[j];
  92.         }
  93. }
  94.     
  95.  
  96. SetPoint(ppa)
  97. PARTICLE *ppa;
  98. {
  99.     ppa->pt.h = xzero + ppa->x[0] / scale;
  100.     ppa->pt.v = yzero - ppa->x[1] / scale;
  101. }
  102.  
  103. Interact(ppa1, ppa2)
  104. PARTICLE *ppa1, *ppa2;
  105. {
  106.     float ftot, f[NDIMS];
  107.     double r2, r;
  108.     double d[NDIMS];
  109.     int i;
  110.  
  111.     r2 = 0;
  112.     for (i = 0 ; i < NDIMS ; i++) {
  113.         d[i] = ppa1->xt[i] - ppa2->xt[i];
  114.         r2 += d[i]*d[i];
  115.     }
  116.     
  117.     r = sqrt(r2);
  118.     if (rmin > r)
  119.         rmin = r;
  120.     ftot = G * ppa1->m * ppa2->m / r2;
  121.     for (i = 0 ; i < NDIMS ; i++) {
  122.         f[i] = ftot * d[i]/r;
  123.         ppa1->f[i] -= f[i];
  124.         ppa2->f[i] += f[i];
  125.     }
  126. }
  127.  
  128. Propagate(ppa, odt, dt)
  129. PARTICLE *ppa;
  130. double odt, dt;
  131. {
  132.     double v, a;
  133.     double step, sesquistep, step2, sesquistep2;
  134.     double v2;
  135.     int i;
  136.     
  137.     step = odt; step2 = 0.5 * odt *odt;
  138.     sesquistep = odt + 0.5 * dt; sesquistep2 = 0.5*sesquistep*sesquistep;
  139.     
  140.     v2 = 0;
  141.     for (i = 0 ; i < NDIMS ; i++) {
  142.         v = ppa->p[i] / ppa->m;
  143.         a = ppa->f[i] / ppa->m;
  144.         ppa->xt[i] = ppa->x[i] + (sesquistep * v) + sesquistep2 * a;
  145.         ppa->x[i] += (step * v) + step2 * a;
  146.         ppa->p[i] += step * ppa->f[i];
  147.         ppa->f[i] = 0;
  148.         v2 += v*v;
  149.     }
  150.  
  151.     if (v2max < v2)
  152.         v2max = v2;
  153. }
  154.  
  155. DrawPath(ppa)
  156. PARTICLE *ppa;
  157. {
  158.     switch (trailstyle) {
  159.         case dot:
  160.             PenMode(erasecode);
  161.             MoveTo(ppa->pt.h, ppa->pt.v);
  162.             LineTo(ppa->pt.h, ppa->pt.v);    /* Keep going */
  163.         case dottrain:
  164.             PenMode(drawcode);
  165.             SetPoint(ppa);
  166.             MoveTo(ppa->pt.h, ppa->pt.v);
  167.             LineTo(ppa->pt.h, ppa->pt.v);
  168.             break;
  169.         case linetrail:
  170.             MoveTo(ppa->pt.h, ppa->pt.v);
  171.             PenMode(drawcode);
  172.             SetPoint(ppa);
  173.             LineTo(ppa->pt.h, ppa->pt.v);
  174.             break;
  175.         default:
  176.             break;
  177.     }
  178. }
  179.  
  180.  
  181. DoIt()
  182. {
  183.     int i, j;
  184.     double t, odt;
  185.     static int tlast;
  186.     int fdraw;
  187.     int whichbutton;
  188.     EventRecord theEvent;
  189.     
  190.     xzero = screenBits.bounds.right/2;
  191.     yzero = screenBits.bounds.bottom/2;
  192.     
  193.     while (TRUE) {
  194.         cpa = 0;
  195.         if (0 == GetParams())
  196.             exit(0);
  197.     
  198.         do {
  199.             whichbutton = GetInit(&rgpa[cpa]);
  200.             if (whichbutton != 3) {            /* if not the END button */
  201.                 cpa++;
  202.             }
  203.         } while (whichbutton == 1);
  204.         if (fcenter)
  205.             CenterParticles(&rgpa, cpa);
  206.         for (i = 0 ; i < cpa ; i++) {
  207.             SetPoint(&rgpa[i]);
  208.         }
  209.         PenNormal();
  210.         PaintRect(&thePort->portRect);
  211.         if (background) {                /* background is not black */
  212.             drawcode = notPatCopy;
  213.             erasecode = patCopy;
  214.         } else {
  215.             InvertRect(&thePort->portRect);
  216.             drawcode = patCopy;
  217.             erasecode = notPatCopy;
  218.         }
  219.         PenMode(drawcode);
  220.         HideCursor();
  221.         
  222.         while (!Button()) {
  223.             if (tlast != (int)(t/drawperiod)) {
  224.                 tlast = t/drawperiod;
  225.                 fdraw = TRUE;
  226.             } else
  227.                 fdraw = FALSE;
  228.             rmin = HUGE;
  229.             for (i = 0 ; i < cpa ; i++) {
  230.                 for (j = i+1 ; j < cpa ; j++) {
  231.                     Interact(&rgpa[i], &rgpa[j]);
  232.                 }
  233.             }
  234.             t += dt;
  235.             odt = dt;
  236.             if (v2max != 0 && rmin != HUGE)
  237.                 dt = rmin/(sqrt(v2max)*precision);
  238.             if (dt > 2*odt)                /* limit growth of timestep    */
  239.                 dt = 2*odt;
  240.             v2max = 0;
  241.             for (i = 0 ; i < cpa ; i++) {
  242.                 Propagate(&rgpa[i], odt, dt);
  243.                 if (fdraw)
  244.                     DrawPath(&rgpa[i]);
  245.             }
  246.             GetNextEvent(0xffff, &theEvent);
  247.         }
  248.         ShowCursor();
  249.     }
  250. }